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Abstract 

It  is  demonstrated  that  finite  precision  Lanczos  and  conjugate  gradient 
computations  for  solving  a  symmetric  positive  definite  linear  system  Ax  = 
b  or  computing  the  eigenvalues  of  A  behave  very  similarly  to  the  exact 
algorithms  applied  to  any  of  a  certain  class  of  larger  matrices.  This  class 
consists  of  matrices  A  which  have  many  eigenvalues  spread  throughout 
tiny  intervals  about  the  eigenvalues  of  A.  The  width  of  these  intervals  is 
a  modest  multiple  of  the  machine  precision  times  the  norm  of  A.  This 
analogy  appears  to  hold,  provided  only  that  the  algorithms  are  not  run 
for  huge  numbers  of  steps.  Numerical  examples  are  given  to  show  that 
many  of  the  phenomena  observed  in  finite  precision  computations  with  A 
can  also  be  observed  in  the  exact  iilgorithms  applied  to  such  a  matrix  A. 


'This  work  was  supported  by  the  Applied  Mathematical  Sciences  Program  of  the  U.S. 
Department  of  Energy  under  contract  DE-AC02-76ER03077. 

'Part  of  this  work  was  performed  while  the  second  author  was  a  visitor  at  the  Courant 
Institute  of  Mathematical  Sciences. 


1      Background  and  Introduction 

The  Lanczos  algorithm  for  computing  eigenvalues  and  eigenvectors  of  a  symmet- 
ric matrix  and  the  conjugate  gradient  algorithm  for  solving  symmetric  positive 
definite  linear  systems  were  introduced  in  the  early  1950's  by  Lanczos  [12]  and 
by  Hestenes  and  Stiefel  [11],  respectively.  It  was  recognized  at  that  time  that 
the  algorithms  often  failed  to  behave  as  they  would  in  exact  arithmetic  due  to 
the  effect  of  rounding  errors.  Engeli,  et.  al.  [5],  for  example,  applied  the  con- 
jugate gradient  method  (without  a  preconditioner)  to  the  biharmonic  equation 
and  demonstrated  that,  for  a  matrix  of  order  n,  convergence  did  not  occur  until 
well  after  step  n  (although  exact  arithmetic  theory  guarantees  that  the  exact 
solution  is  obtained  after  n  steps).  For  this  and  other  reasons,  the  algorithms 
did  not  gain  widespread  popularity  at  that  time. 

With  the  idea  of  preconditioning  in  the  conjugate  gradient  method,  interest 
in  this  algorithm  was  revived  in  the  early  1970's,  with  several  important  papers 
appearing,  including  work  by  Reid  [16]  and  by  Concus,  Golub,  and  O'Leary  [4]. 
News  of  the  effectiveness  of  the  conjugate  gradient  method  as  an  iterative  tech- 
nique quickly  spread  throughout  the  scientific  computing  community,  and  the 
algorithm  soon  became  the  most  popular  method  for  solving  symmetric  positive 
definite  linear  systems.  Although  the  effect  of  rounding  errors  on  the  conjugate 
gradient  algorithm  was  not  well  understood,  it  was  observed  numerically  that 
(with  a  good  preconditioner)  either  the  method  converged  before  rounding  er- 
rors had  any  significant  effect  on  the  iterates,  or,  whatever  the  effect  of  roundofi", 
it  was  not  catastrophic. 

Further  attempts  were  made  to  understand  the  effect  of  rounding  errors  on 
these  two  mathematically  equivalent  algorithms,  and  why,  if  they  were  run  for 
enough  steps  to  be  significantly  affected  by  roundoff,  the  effects  were  not  disas- 
trous. Wozniakowski  [19]  considered  a  special  version  of  the  conjugate  gradient 
(CG)  algorithm  and  showed,  essentially,  that  a  finite  precision  implementation 
converged  at  least  as  rapidly  as  the  method  of  steepest  descent.  More  precisely, 
if  the  linear  system  to  be  solved  is  Ax  =  6,  if  e*  =  x  —  x*  denotes  the  error  in 
the  /b"*  iterate,  and  if  k  is  the  condition  number  of  A,  then  the  >4-norm  of  the 
error  at  step  k  satisfies 

\\e'\\A    <   (l+20(^)||e*-MU  +  0(0,  (1) 

where  e  is  the  unit  roundoff  of  the  machine  and  0(f)  denotes  terms  involving 
the  product  of  i  with  the  norm  of  various  powers  of  A,  ar*,  and  a  constant. 
Wozniakowski  also  gave  a  bound  on  the  ultimately  attainable  accuracy  using 
this  special  version  of  the  CG  algorithm: 

limsup  11411^    <   ,^1'C.  (2) 


Iimsupll^;i4^    <e.\\A\\C. 

l^oo  11^    II2 

Cullum  and  Willoughby  [3]  proved  a  result  similar  to  (1)  for  a  more  standard 
version  of  the  CG  algorithm. 

The  bound  (1)  is  a  large  overestimate  of  the  actual  error,  however.  If  the 
CG  algorithm  really  converged  as  slowly  as  the  method  of  steepest  descent,  it 
would  seldom  be  used.  Methods  such  as  the  Chebyshev  algorithm  or  Richard- 
son's method  would  be  far  superior.  The  bound  (1)  was  derived  by  considering 
individual  steps  of  the  CG  algorithm  -  assuming  only  that  a  particular  step  k 
is  implemented  accurately,  it  follows  that  the  error  at  step  k  is  reduced  at  least 
as  much  as  it  would  be  by  a  steepest  descent  step.  Yet,  an  example  due  to 
Crowder  and  Wolfe  [2]  shows  that  unless  one  considers  all  steps  of  the  CG  al- 
gorithm, one  cannot  hope  to  establish  much  faster  convergence  than  this.  That 
is,  if  the  initial  search  direction  is  chosen  incorrectly  but  all  other  steps  of  the 
algorithm  are  implemented  exactly,  then  convergence  may  be  almost  as  slow  as 
the  method  of  steepest  descent. 

The  following  simple  3  by  3  example  presented  in  [2]  illustrates  this  phe- 
nomenon. 


V6\      0       r  4v/30  \        3^6 


If  r°  is  the  initial  residual  for  a  linear  system  with  coefficient  matrix  A,  and  if, 
instead  of  taking  the  initial  search  direction  p°  to  be  r°  we  set  p°  as  above,  then 
if  the  remaining  CG  formulas, 

'  *  ^        p'-\      it  =  1,2,..., 


^     ~  p^-i^>lp*-' 

are  implemented  exactly,  then  r*  will  satisfy 


,  /  1         0  0        \ 

r*    =    -Qr''-\       it  =  1,2,...,       where  Q   =         0      -1/5      -2n/6/5       , 
^  V  0     2v/6/5       -1/5     / 

and  the  j4-norm  of  the  error  will  be  reduced  by  a  factor  of  |  at  each  iteration. 
This  is  somewhat  faster  than  the  steepest  descent  bound. 


K-  1    _    _9_ 
K+l    ~    IT' 


but  it  is  slower  than  the  Chebyshev  method  which  would  converge  at  an  asymp- 
totic rate 

#4   «    .52. 

While  most  work  on  the  CG  algorithm  was  focusing  on  individual  steps,  the 
effects  of  rounding  errors  on  the  Lanczos  algorithm  were  being  studied  from  a 
more  global  point  of  view  -  considering  the  effects  of  roundoff  over  the  entire 
course  of  the  computation.  Paige  [13]  wrote  the  Lanczos  recurrence  equations 
in  matrix  form  along  with  the  matrix  of  perturbations  resulting  from  finite  pre- 
cision arithmetic.  Using  this  formulation,  he  analyzed  the  loss  of  orthogonality 
among  the  Lanczos  vectors.  He  showed  that  loss  of  orthogonality  is  only  in 
the  direction  of  converged  Ritz  vectors.  From  this  it  followed  that  at  \easi  one 
eigenpair  must  converge  by  step  n.  Paige  later  showed  also  that  Ritz  values 
"stabilize"  only  to  points  near  eigenvalues  of  A  [14].  The  implications  of  these 
results  as  far  as  the  rate  of  convergence  of  the  Lanczos  or  CG  algorithm  were 
not  so  clear. 

Parlett  and  Scott  [15]  used  Paige's  results  to  suggest  a  "fix"  for  the  Lanczos 
algorithm  -  selective  orthogonalization.  This  requires  saving  the  Lanczos  vec- 
tors and  orthogonalizing  against  Ritz  vectors  as  they  converge.  Further  work 
on  reorthogonalization  strategies,  as  well  as  on  understanding  the  behavior  of 
the  algorithms  without  reorthogonalization,  was  carried  out  by  Simon  [17].  He 
showed  that  until  approximate  orthogonality  is  lost,  the  tridiagonal  matrix  gen- 
erated by  a  finite  precision  Lanczos  computation  is,  indeed,  the  approximate 
projection  of  the  matrix  A  onto  the  span  of  the  Lanczos  vectors  (which  may 
or  may  not  be  the  desired  Krylov  space).  Grcar  [8]  attempted  a  forward  error 
analysis  of  the  conjugate  gradient  algorithm.  He  showed  that  under  a  certain  as- 
sumption, called  the  "projection  property",  the  coefficients  generated  in  a  finite 
precision  CG  computation  are  within  about  e  of  those  that  would  be  generated 
by  the  exact  algorithm,  as  long  as  the  vectors  remain  within  about  ^/e  of  the 
exact  ones.  Thus  the  initial  deviation  from  exact  arithmetic  can  be  analyzed  as 
if  the  coefficients  were  given  rather  than  computed  at  each  step. 

In  [10]  a  form  of  backward  error  analysis  was  developed  for  the  Lanczos  and 
CG  algorithms.  There  it  was  shown  that  finite  precision  computations,  run  for 
no  more  than  some  number  J  steps,  generate  the  same  tridiagonal  matrices  at 
each  step  as  the  exact  algorithms  applied  to  a  larger  matrix  A,  having  possibly 
many  more  eigenvalues  than  A,  but  whose  eigenvalues  all  lie  within  tiny  intervals 
about  the  eigenvalues  of  A.  A  bound  on  the  size  of  these  intervals  was  derived  in 
terms  of  the  machine  precision  e  and  the  bound  J  on  the  number  of  steps.  Using 
this  analogy  it  was  possible,  in  some  cases,  to  derive  more  interesting  bounds  on 
the  convergence  rate  of  the  CG  and  Lanczos  algorithms.  For  example,  assuming 
that  the  widths  of  the  intervals  containing  the  eigenvalues  of  A  are  much  smaller 
than  the  smallest  eigenvalue  of  A,  it  follows  that  the  condition  number  of  j4  will 
be  approximately  the  same  as  that  of  A.  Hence  any  exact  arithmetic  error  bound 


in  terms  of  the  condition  number  k  o{  A  will  hold,  to  a  close  approximation,  for 
finite  precision  computations;  e.g.. 

This  error  bound  was  conjectured  in  [19],  but  it  could  not  be  proved  with  the 
approach  used  there.  The  bound  derived  in  [10]  on  the  size  of  these  intervals 
was,  unfortunately,  a  large  overestimate.  Thus,  while  bounds  such  as  (3)  were 
established  if  the  machine  precision  (  was  small  enough,  in  many  realistic  prob- 
lems the  bound  on  the  interval  size  was  too  large  to  provide  useful  information. 
A  procedure  was  given  to  actually  compute  the  matrix  A,  however,  and  numeri- 
cal computation  of  this  matrix  for  many  examples  indicates  that  the  eigenvalues 
of  A  are  actually  contained  in  much  smaller  intervals  than  the  proven  bound 
would  suggest. 

In  this  paper  we  demonstrate  numerically  that  the  behavior  of  finite  preci- 
sion Lanczos  and  CG  computations  is  not  only  identical  to  that  of  the  exact 
algorithms  applied  to  a  particular  matrix  A,  but  also  very  similar  to  that  of 
the  exact  algorithms  applied  to  any  matrix,  say.  A,  which  has  many  eigenvalues 
spread  throughout  tiny  intervals  about  the  eigenvalues  of  A.  The  size  of  the 
intervals  is  a  modest  multiple  of  the  machine  precision.  It  is  not  clear  if  this 
similarity  is  maintained  if  the  algorithms  are  run  for  huge  numbers  of  steps  (say, 
^  10^),  but  for  more  realistic  computations,  the  similarity  is  demonstrated.  For 
test  problems,  we  consider  a  class  of  matrices  introduced  in  [18].  There  the  be- 
havior of  the  finite  precision  computations  was  compared  with  exact  arithmetic 
theory  and  shown  to  give  surprising  results.  In  this  paper  we  show  why  this 
behavior  is  to  be  expected. 

Finite  precision  CG  computations  for  solving  an  n  by  n  symmetric  positive 
definite  linear  system  Ax  =  b  sometimes  fail  to  converge  after  n  steps,  especially 
when  n  is  small.  In  such  caises,  it  is  demonstrated  that  exact  CG  applied  to  the 
corresponding  large  linear  system  Ax  =  6  also  requires  more  than  n  iterations 
to  converge.  More  commonly,  finite  precision  CG  computations  converge  in  far 
fewer  than  n  steps,  and  the  same  holds  for  the  exact  CG  algorithm  applied  to  any 
matrix  A  whose  eigenvalues  are  clustered  in  tiny  intervals  about  the  eigenvalues 
of  A.  Frequently,  finite  precision  CG  computations  go  through  several  steps  at 
which  there  is  only  a  modest  reduction  in  the  error  and  then  at  the  next  step 
there  is  a  very  sharp  decrease  in  the  error.  This  same  behavior  is  observed  in 
the  exact  CG  algorithm  applied  to  matrices  A  whose  eigenvalues  are  distributed 
in  n  tight  clusters  about  the  eigenvalues  of  A. 

Related  to  this  phenomenon  of  slow  convergence  followed  by  a  sudden  drop 
in  the  error,  is  the  phenomenon  of  multiple  "copies"  of  eigenvalues  appearing  in 
finite  precision  Lanczos  computations.  Finite  precision  Lanczos  computations 
frequently  generate  several  close  approximations  to  some  of  the  eigenvalues  of  A 
before  finding  any  close  approximations  to  some  of  the  other  eigenvalues.  Anal- 


ogously,  depending  on  how  the  clusters  of  the  larger  matrix  A  are  distributed, 
the  exact  Lanczos  algorithm  applied  to  A  may  find  several  eigenvalues  with- 
in some  of  the  clusters  before  finding  any  in  some  of  the  other  clusters.  It  is 
demonstrated  that  the  rate  of  occurrence  of  multiple  "copies"  of  eigenvalues  in 
finite  precision  Lanczos  computations  with  matrix  A  is  very  similar  to  the  rate 
at  which  the  exact  Lanczos  algorithm  applied  to  A  finds  different  eigenvalues 
within  the  same  cluster. 

The  following  sections  2-4  present  numerical  examples  to  demonstrate  these 
phenomena.  Implications  of  this  analogy  are  discussed  in  section  5. 

2      Description  of  Numerical  Experiments 

The  matrices  considered  in  this  paper  were  introduced  in  [18]  and  have  eigen- 
values of  the  form 

A.  =  A, +  -?— ^(A„-AiK-*,       f=2,...,n,       pe(0,l),  (4) 

n  —  1 

where  n,  Ai,  and  k  -  A„/Ai  are  fixed.  For  most  of  our  experiments  we  have 
taken  n  =  24,  Ai  =  .1,  /c  =  1000,  and  p=  .4,  .6,  .8,  .9,  1.0.  These  eigenvalues 
are  plotted  for  each  value  of  p  in  Fig.  1.  Eigenvalues  that  are  too  close  to 
be  distinguished  on  the  horizontal  axis  have  been  plotted  vertically.  For  the 
smaller  />-values,  the  eigenvalues  are  very  tightly  clustered  at  the  lower  end  of 
the  spectrum.  The  minimal  difference,  Aj  -  Ai,  corresponding  to  p  =  .4,  .6,  .8, 
.9,  and  1.0  is  7.6e  -  9,  5.7e  -  5,  .032,  .43,  and  4.3,  respectively. 

The  algorithm  used  in  these  experiments  for  solving  a  symmetric  positive 
definite  linear  system  Ax  =  b  and  computing  the  eigensystem  of  v4  is  as  follows 
[11,12]: 

Given  an  initial  guess  x°,  compute  r°  =  6  -  Ax° ,  and  set  p°  =  r°. 
For  it  =  1,2,... 


Compute  aic-i  = 


p- 


Take  i*  =  x*"^  -|-ai_ip*"^ 
Compute  r*  =  r*"'  -  at^iAp^'^. 

Compute  0k  -  ^tliT^t,.- 

Set  Tt.i+i  =  Ti+i,i  =  ^^. 
Takep*  =  r*4-/3tp*-'. 

The  tridiagonal  matrix  T  generated  at  step  it  will  be  denoted  T^  \  and  its  eigen- 
values are  taken  as  approximate  eigenvalues  of  A  (Ritz  values).  The  eigenvectors 
o{  A  can  also  be  approximated  if  the  previous  residual  vectors,  r°,...,r*~',  have 


been  saved,  but  we  will  not  discuss  the  computation  of  eigenvectors  here.  When 
solving  linear  systems,  we  will  refer  to  this  algorithm  as  the  conjugate  gradient 
method,  or,  CG,  while  when  using  it  to  compute  eigenvalues  we  will  refer  to  it 
as  the  Lanczos  algorithm.  The  equivalence  of  this  method  to  the  usual  Lanczos 
process  [12]  in  exact  arithmetic  is  well-known,  and  arguments  in  [3,10]  establish 
similar  behavior  in  finite  precision  arithmetic  as  well.  Numerical  experiments 
with  other  variants  of  this  algorithm  have  yielded  similar  results,  as  described 
in  [18]. 

The  above  algorithm  was  applied  to  matrices  A  of  the  form 

A  =  UAU'^, 

where  t/  is  a  random  orthogonal  matrix  and  A  =  diag{\i, ...,  An)  is  defined  in 
(4).  In  all  cases,  a  random  right-hand-side  vector  and  a  zero  initial  guess  were 
used.  Experiments  were  carried  out  on  a  Sun  Sparcstation  using  double  precision 
arithmetic  (about  16  decimal  digits).  Most  of  the  experiments  were  performed 
using  MATLAB,  making  the  algorithms  relatively  simple  to  implement. 

We  first  compare  finite  precision  computations  for  solving  Ax  =  6  or  com- 
puting the  eigenvalues  of  A  to  the  exact  arithmetic  algorithms  applied  to  the 
same  problem.  "Exact  arithmetic"  was  simulated  by  saving  all  residual  vectors 
and  using  full  reorthogonalization  at  each  step.  That  is,  the  formula  for  r* 
becomes 


For  kount  - 

■■  1,2, 

iV 

For  j  = 

l,it- 

-  1, 

Endfor 

-.T^- 

77^ 

Endfor 

It  is  shown  in  [14]  that  the  iterates  generated  using  this  modified  algorithm  do, 
indeed,  resemble  those  that  would  be  generated  by  the  exact  algorithm  applied 
to  a  slightly  different  matrix  (of  the  same  order)  with  a  slightly  diflTerent  initial 
vector. 

We  also  compare  the  finite  precision  computations  involving  the  matrix  A 
to  "exact  arithmetic"  (full  reorthogonalization)  computations  involving  a  larger 
matrix  j4.  The  matrix  A  was  taken  to  have  a  total  of  1  In  eigenvalues,  with  eleven 
eigenvalues  uniformly  distributed  throughout  each  of  n  tiny  intervals  about  the 
eigenvalues  of  j4.  Several  of  the  experiments  were  also  performed  with  a  matrix 
A  having  twenty  one  eigenvalues  evenly  distributed  in  each  of  these  same  inter- 
vals, and  the  results  were  identical,  to  within  graphical  accuracy.  Most  of  the 
experiments  were  performed  with  intervals  of  width  10~^^,  or,  approximately 

bOc\\A\\ 


where  e  =  2~^^  is  the  unit  roundoff  of  the  machine.  Some  experiments  were 
performed  with  different  size  intervals  to  see  the  effect  of  the  interval  width  on 
the  behavior  of  the  algorithms.  The  finite  precision  computation  for  Ax  =  6, 
with  initial  guess  x°,  or,  initial  residual  r°,  was  compared  to  the  exact  arithmetic 
computation  for  Ax  =  6,  with  initial  guess  i",  initial  residual  f°,  where 

A  =  diag{Xi^i,...,Xi,m,  A2,i,...,A2,m,  ■••,  -^n.i,  ■  •-,  A^  „,),  (5) 


m+l 


Xii=Xi  +  - ^—-6,       j=l,...,m,       m=ll,      6  =  10"^- 

m  —  1 

f°  =  b,      6,,i  =  6.,2=...  =  6,,m.    and    Yl^b.^jf  ^  (U'^bf-,       i=l,...,n, 

;  =  i 

where  Ai An  are  the  eigenvalues  of  .4  and  U  is  the  orthonormal  matrix  of 

eigenvectors  of  A. 

3      Results  of  CG  Computations 

Here  we  show  the  remarkable  similaritj'  between  a  finite  precision  CG  compu- 
tation to  solve  Ax  =  b,  with  initial  residual  r° ,  and  and  the  exact  CG  algorithm 
applied  to  the  larger  linear  system  Ax  =  b,  with  initial  residual  r°. 

Fig.  2a  shows  the  convergence  of  finite  precision  CG  computations  applied 
to  the  linear  system  Ax  —  b,  for  the  five  different  p- values  listed  in  the  previous 
section.  The  right-hand  side  vector  6  was  taken  to  have  random  components, 
uniformly  distributed  between  0  and  1,  and  the  initial  guess  i°  was  taken  to  be 
zero.  The  .4-norm  of  the  error  at  each  iteration  divided  by  the  v4-norm  of  the 
initial  error, 

<x-x^,A{x-x^)>^f^ 
<x-x°,A{x-x°)  >i/2 

is  plotted.  (The  "exact"  solution  x  was  computed  as  UA~^U^b.)  Note  that 
although  exact  arithmetic  theory  ensures  that  the  correct  solution  is  obtained 
after  n  =  24  steps,  the  finite  precision  computation  requires  significantly  more 
than  n  steps  for  some  of  the  />-values.  For  certain  values  of  p  the  computation 
seems  to  be  considerably  more  affected  by  rounding  errors  than  for  other  values. 
Convergence  slows  as  p  goes  from  .4  to  .6  to  .8,  but  then  improves  as  p  reaches 
.9  and  1. 

For  comparison,  Fig.  2b  shows  the  convergence  of  the  exact  CG  algorithm 
applied  to  the  same  linear  system,  with  the  same  initial  guess.  In  contrast  to 
the  finite  precision  computation,  there  is  little  difference  between  the  results  for 
p  =  .8,  .9,  and  1.0.  with  the  slowest  exact  arithmetic  convergence  rate  occurring 
for  p  =  .9. 


The  convergence  of  the  finite  precision  computations  much  more  closely  re- 
sembles that  of  the  exact  CG  algorithm  applied  to  the  linear  system  Ax  =  6 
(defined  in  (5)),  shown  in  Fig.  2c.  Here  the  j4-norm  of  the  error  at  each  iter- 
ation divided  by  the  A-norm  of  the  initial  error  is  plotted.  As  with  the  finite 
precision  CG  computations,  convergence  slows  as  p  goes  from  .4  to  .6  to  .8,  but 
then  improves  for  p  =  .9  and  1.  The  qualitative  convergence  behavior  in  Fig. 
2c  is  also  similar  to  that  of  the  finite  precision  CG  computations,  in  that,  for 
certain  p  values,  both  go  through  stages  at  which  little  improvement  is  made 
for  several  steps  and  then  a  sharp  drop  in  the  error  is  seen  at  a  subsequent  step. 

To  see  the  effect  of  the  interval  size  on  the  convergence  rate  of  the  exact  CG 
algorithm  applied  to  a  matrix  A  with  eigenvalues  clustered  in  these  intervals,  we 
tried  several  different  interval  sizes  for  the  case  p  =  .6.  That  is,  we  considered 
matrices  A  whose  eigenvalues  were  clustered  in  intervals  of  width  6  =  10"^^, 
10~'^,  10~^°,  and  10~^  about  the  eigenvalues  of  A.  As  before,  each  matrix  A 
was  taken  to  have  eleven  eigenvalues  in  each  interval,  uniformly  distributed,  and 
the  initial  residuals  r°  were  set  according  to  (5).  Fig.  3  shows  the  convergence 
of  the  exact  CG  algorithm  applied  to  the  different  problems  Ax  =  6,  along 
with  that  of  the  finite  precision  computation  for  Ax  =  6.  While  the  exact 
computation  with  the  matrix  of  interval  width  10"^^  most  closely  resembles  the 
finite  precision  computation,  similar  qualitative  behavior  is  reflected  in  all  of 
these  computations,  except  perhaps  the  one  with  interval  width  10~^,  which 
is  considerably  slower  to  converge.  Thus,  it  appears  that  a  precise  estimate 
of  this  interval  width  is  not  even  necessary  to  predict  the  qualitative  behavior 
of  finite  precision  CG  computations.  They  resemble  exact  CG  computations 
for  any  matrix  A  with  eigenvalues  spread  throughout  small  intervals  about  the 
eigenvalues  of  A,  and  the  interval  size  can  be  anywhere  within  a  rather  wide 
range.  The  remainder  of  the  CG  comparisons  will  use  10~^^  as  the  interval 
width  for  the  eigenvalues  of  A. 

While  most  of  our  experiments  have  been  performed  with  very  small  matrices 
(n  =  24),  similar  phenomena  can  be  observed  with  larger  matrices,  for  which 
the  CG  algorithm  is  more  often  used.  For  large  matrices,  our  comparisons  with 
exact  arithmetic  computations  for  A  become  quite  time-  and  storage-consuming, 
however,  since  A  is  eleven  times  as  large  as  A  and  it  is  necessary  to  save  all 
residual  vectors  and  reorthogonalize  at  every  step.  Still,  we  have  performed  one 
experiment  with  a  matrix  A  of  order  100.  The  eigenvalues  of  ^  are  still  defined 
by  formula  (4),  with  n  =  100,  Aj  =  .1,  k  =  1000,  and  we  took  p  =  .8.  Fig.  4 
shows  the  convergence  of  finite  precision  CG  for  solving  Ax  =  6,  exact  CG  for 
solving  Ax  =  6,  and  exact  CG  for  solving  the  larger  problem  Ax  =  6.  Note  the 
close  resemblance  between  the  first  and  last  of  these  curves  and  the  significant 
differences  from  the  exact  CG  computation  for  Ax  =  6.  Although  this  type  of 
behavior  is  frequently  seen  in  practice,  it  may  not  be  realized  that  rounding 
errors  are  significantly  affecting  the  convergence  rate,  since  there  is  no  exact 
arithmetic  computation  with  which  to  compare. 


4     Results  of  Lanczos  Computations 

Here  we  show  the  similarity  between  the  eigenvalue  approximations  generated 
at  each  step  of  a  finite  precision  Lanczos  computation  with  matrix  A  and  initial 
vector  r°  and  those  generated  at  each  step  of  an  exact  Lanczos  computation 
with  the  larger  matrix  A  and  initial  vector  r°,  defined  in  (5). 

Figs.  5-7a  show  the  eigenvalue  approximations  generated  by  a  finite  precision 
Lanczos  computation  with  matrix  A,  for  the  cases  p  =  .6,  .8,  1.0.  Similar 
results  were  observed  with  the  other  p-values.  In  order  to  distinguish  clustered 
eigenvalues,  we  have  plotted  on  the  vertical  Eixis  not  the  actual  eigenvalues,  but 
the  index  of  each  eigenvalue  of  A,  from  1  to  n.  An  approximate  eigenvalue  that 
lies,  say,  1/r  of  the  way  between  eigenvalue  i  and  eigenvalue  j '+  1  of  A,  will  be 
plotted  on  the  graph  at  y- value  i+  ^-  Eigenvalue  approximations  that  are  too 
close  to  be  distinguished  on  the  vertical  axis  have  been  plotted  horizontally.  For 
clarity,  we  have  plotted  the  eigenvalue  approximations  only  at  every  fourth  step. 
In  most  cases,  we  see  multiple  copies  of  the  larger  eigenvalues  appearing  before 
any  close  approximations  to  the  smaller,  clustered  eigenvalues  appear.  It  should 
be  remembered,  however,  that  for  the  smaller  p-values,  these  small  eigenvalues 
are  very  tightly  clustered,  and  there  may  be  a  close  approximation  to  the  cluster, 
even  though  the  individual  eigenvalues  have  not  been  identified.  Only  in  the 
case  p  =  1  does  the  finite  precision  computation  find  all  n  eigenvalues  by  step 
n. 

For  comparison.  Figs.  5-7b.  show  the  eigenvalue  approximations  generated 
every  fourth  step  of  the  exact  Lanczos  algorithm  applied  to  the  same  matrices, 
with  the  same  initial  vectors.  Here  no  "multiple  copies"  are  observed,  and  all 
of  the  eigenvalues  are  identified  after  n  steps. 

The  eigenvalue  approximations  generated  by  the  finite  precision  computa- 
tions much  more  closely  resemble  those  generated  by  the  exact  Lanczos  algcn 
rithm  applied  to  the  matrix  A,  with  initial  vector  r°,  shown  in  Figs.  5-7c.  In 
these  figures  we  again  see  multiple  close  approximations  to  the  larger  eigenval- 
ues appearing  before  step  n,  except  in  the  case  p  =  1.  The  rate  of  appearance 
of  these  multiple  copies  also  appears  to  be  similar  to  that  in  the  finite  precision 
computations  with  matrix  A. 

We  point  out  that  for  p  =  1.0,  the  effect  of  roundoff  on  the  Lanczos  process 
is  minimal  (cf.  Fig.  7a,b;  also  Fig.  2a,b).  After  n  iterations  the  process  is 
"restarted"  and  it  computes  all  eigenvalues  twice  in  2n  iterations.  For  p  < 
1.0,  the  "restarting"  is  more  frequent  and  multiple  copies  of  large  eigenvalues 
are  computed  simultaneously  with  single  approximations  to  small  eigenvalues. 
It  can  be  observed  that  if  a  finite  precision  Lanczos  computation  generates  a 
close  approximation  to  each  eigenvalue  of  A  at  some  step,  then  the  error  in  the 
corresponding  finite  precision  CG  computation  drops  dramatically  at  that  step. 
See  Fig.  6a  and  Fig.  2a  (p  =  .8),  between  32  and  36  iteration  steps. 

Again,  to  see  the  effect  of  the  interval  width  on  the  eigenvalue  approximation- 
s  generated  by  an  exact  Lanczos  computation  for  a  matrix  A  with  eigenvalues 
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clustered  in  these  intervals,  we  tried  several  diflferent  interval  sizes,  for  the  case 
p  =  -6.  That  is,  we  considered  matrices  A  whose  eigenvalues  were  clustered 
in  intervals  of  width  6  =  10~'^,  10~*°,  and  10"^  about  the  eigenvalues  of  A. 
As  before,  each  matrix  A  was  taken  to  have  eleven  eigenvalues  in  each  inter- 
val, uniformly  distributed,  and  the  initial  residuals  f°  were  set  according  to  (5). 
Figs.  8a-c  show  the  eigenvalue  approximations  generated  at  every  fourth  step 
by  the  exact  Lanczos  algorithm  applied  to  each  of  these  matrices  A.  Note  the 
similarity  between  each  of  these  figures  (as  well  as  Fig.  5c  with  6  =  10"^^)  and 
Fig.  5a,  showing  the  eigenvalue  approximations  generated  by  a  finite  precision 
Lanczos  computation  for  the  matrix  A.  Fig.  8a  (6  =  10"^^)  shows  the  closest 
resemblance  to  the  finite  precision  computation,  while  Fig.  8c  [6  =  10"^)  has 
somewhat  more  copies  of  the  larger  eigenvalues  and  somewhat  fewer  approxi- 
mations to  the  smaller  ones. 

Finally,  in  Figs.  9a-c,  we  have  plotted  results  from  a  larger  problem,  with  n  = 
100,  p  =  .8.  Again  note  the  similarities  between  the  eigenvalue  approximations 
generated  by  the  finite  precision  Lanczos  computation  for  the  matrix  A  (Fig. 
9a)  and  those  generated  by  the  exact  arithmetic  Lanczos  computation  for  the 
matrix  A  (Fig.  9c).  Unlike  the  exact  Lanczos  computation  for  A  (Fig.  9b),  these 
procedures  both  generate  multiple  close  approximations  to  some  of  the  larger 
eigenvalues  before  finding  any  close  approximations  to  some  of  the  smaller  ones. 
The  rate  at  which  these  multiple  approximations  appear  is  also  similar  in  Figs. 
9a  and  c. 


5      Further  Discussion  and  Open  Questions 

We  have  demonstrated  numerically  that  the  behavior  of  finite  precision  Lanczos 
and  CG  computations  with  a  matrix  A  closely  resembles  that  of  the  exact  algo- 
rithms applied  to  matrices  A  which  have  many  eigenvalues  spread  throughout 
tiny  intervals  about  the  eigenvalues  o{  A.  In  [10]  it  was  proved  that  the  eigen- 
value approximations  generated  by  a  finite  precision  Lanczos  computation  are 
identical  io  those  generated  by  the  exact  algorithm  applied  to  a  certain  larger 
matrix  A,  and  that  the  ^-norm  of  the  error  in  the  equivalent  finite  precision  CG 
computation  is  reduced  at  approximately  the  same  rate  as  the  ^-norm  of  the 
error  in  the  exact  algorithm.  Eigenvalues  of  the  matrix  A  lie  in  tiny  intervals 
about  the  eigenvalues  of  A  (provided  that  the  finite  precision  computation  is 
not  run  for  too  many  steps),  but  A  might  have  many  or  only  a  few  eigenvalues 
in  some  of  these  intervals. 

Using  these  analogies,  the  problem  of  estimating  and  bounding  the  conver- 
gence rates  of  these  algorithms  in  finite  precision  arithmetic  reduces  to  a  problem 
of  estimating  or  bounding  the  convergence  rates  of  the  exact  algorithms  applied 
to  certain  classes  of  matrices.  If  an  error  bound  can  be  established  for  the  exact 
algorithm  applied  to  every  matrix  whose  eigenvalues  lie  within  these  interval- 
s,  then  it  will  hold  (at  least  to  a  close  approximation)  for  the  finite  precision 
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computation.  If  an  error  estimate  is  good  for  the  exact  algorithms  applied  to 
every  matrix  with  many  eigenvalues  spread  throughout  these  intervals,  then  it 
will  also  be  a  good  estimate  for  the  finite  precision  computation.  We  will  not 
derive  such  bounds  and  estimates  here,  but  it  is  not  difficult  to  see  how  they 
might  be  derived  (as,  for  example,  in  (3)),  and  some  examples  are  given  in  [10]. 
A  question  that  is  frequently  asked  is  whether  a  finite  precision  Lanczos 
computation  eventually  finds  all  eigenvalues  of  a  matrix.  Using  the  analogy 
developed  in  [10]  between  finite  precision  Lanczos  computations  run  for  no  more 
than  J  steps  and  the  exact  algorithm  applied  to  a  matrix  whose  eigenvalues  lie 
within  intervals  of  width,  say,  6j ,  about  the  eigenvalues  of  A,  this  question  can 
be  translated  as  follows:  Is  there  a  J  such  that  the  exact  Lanczos  algorithm 
applied  to  every  matrix  A  whose  eigenvalues  lie  within  intervals  of  width  6j 
about  the  eigenvalues  of  >4  -  with  initial  vector  r°  satisfying 

X:(r°,u''V    «    (r^«')^       f=l,...,n,  (6) 

I 
where  u\  ...,u"  are  the  eigenvectors  of  yl  and  u''',  (.  =  1, ...,  are  the  eigenvectors 
of  j4  corresponding  to  the  eigenvalues  clustered  about  A,-,  »  =  l,...,n  -  finds  at 
least  one  eigenvalue  from  each  cluster  within  J  steps?  We  know  of  no  simple 
and  general  sufficient  conditions  for  the  existence  of  such  a  number  J ,  so  this 
question  remains  open. 

Of  course,  if  the  interval  widths  6j  could  be  bounded  (with  a  suitably  small 
bound)  independent  of  J,  then  it  would  follow  that  a  finite  precision  Lanc- 
zos computation  would  eventually  find  every  eigenvalue  of  A  whose  eigenvector 
contained  a  nonnegligible  component  in  the  initial  vector.  For  (6)  implies,  in 
this  case,  that  the  interval  about  each  eigenvalue  has  a  nonzero  weight.  ^From 
Favard's  theorem  [6]  it  would  follow  that  the  characteristic  polynomials  of  the 
tridiagonal  matrices  generated  by  a  finite  precision  Lanczos  computation  were 
the  orthogonal  polynomials  for  a  certain  measure  whose  support  lies  in  the  union 
of  intervals  (J"_][A,  —  (5,  A,  +  6],  where  b  is  the  bound  on  the  interval  size.  But 
the  roots  of  the  orthogonal  polynomials  are  known  to  converge  to  all  weighted 
points  (see,  for  example  [1]),  so  they  would  converge  to  at  least  one  point  in 
each  of  these  intervals.  Thus,  we  could  then  conclude  that  a  finite  precision 
Lanczos  computation  would  eventually  find  an  approximation  within  6  of  each 
eigenvalue  of  A. 

Whether  this  interval  size  can  be  bounded  by  a  small  number  8  independent 
of  J  is  an  open  question.  (Of  course,  there  is  some  bound  that  is  independent 
of  J.  From  Gershgorin's  theorem  and  the  formulas  for  the  elements  of  the  tridi- 
agonal matrices,  it  can  be  seen  that  all  approximate  eigenvalues  generated  by  a 
finite  precision  Lanczos  computation  lie  in  the  interval  [\min  —  "^^maz ,^^maz\- 
Hence  the  measure  for  which  the  characteristic  polynomials  are  orthogonal  has 
its  support  in  this  set.  But  this  is  not  a  very  interesting  bound.) 

To  try  and  determine  whether  such  a  bound  8  exists,  we  have  run  a  finite 
precision  Lanczos  computation  for  the  case  p  =  .6,  n  =  24,  for  10^  steps,  and  we 
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have  computed  the  spread  of  the  eigenvalue  approximations  clustered  about  the 
largest  eigenvalue  of  A.  That  is,  we  have  taken  all  eigenvalue  approximations 
which  are  closer  to  the  largest  eigenvalue  of  A  than  to  any  other  eigenvalue  o(  A, 
and  we  have  computed  the  difference  between  the  largest  of  these  and  the  second 
smallest  of  these.  By  the  interlace  theorem,  every  future  tridiagonal  matrix  will 
have  an  eigenvalue  greater  than  the  largest  of  these  approximations  and  an 
eigenvalue  between  each  pair  of  these  approximations,  and  so,  this  difference 
gives  a  lower  bound  on  the  interval  size  6  in  which  the  weighted  points  lie. 

Results  using  double  and  single  precision  arithmetic  are  plotted  in  Figs.  10a 
and  b.  For  these  experiments,  we  used  the  standard  formulation  of  the  Lanczos 
algorithm,  rather  than  the  CG  form  presented  earlier,  to  avoid  problems  with 
underflow.  As  can  be  seen  from  the  figure,  this  lower  bound  continues  to  grow 
with  J  at  least  out  to  J  =  10^,  but  it  is  growing  very  slowly.  Whether  it  stops 
growing  at  some  value  significantly  less  than  the  distance  to  the  next  largest 
eigenvalue,  or  whether  these  eigenvalue  approximations  would  eventually  fill  the 
entire  Gershgorin  interval,  we  cannot  say.  This  remains  an  open  question. 

6      Conclusions 

We  have  found  the  analogy  between  finite  precision  CG/Lanczos  computations 
and  the  exact  algorithms  applied  to  a  larger  matrix  with  nearby  eigenvalues  to 
be  useful  in  understanding  and  predicting  the  behavior  of  such  computations. 
The  proven  identity  between  finite  precision  computations  for  A  and  exact  com- 
putations for  A  enables  one  to  prove  results  about  finite  precision  CG/Lanczos 
computations.  The  demonstrated  similarity  between  finite  precision  computa- 
tions for  A  and  exact  computations  for  matrices  A  enables  one  to  estimate  the 
actual  behavior  of  finite  precision  CG/Lanczos  computations.  It  provides  a  nice 
explanation  of  the  phenomena  observed  in  [18],  for  example.  The  proven  bound 
on  the  size  of  the  intervals  containing  the  eigenvalues  of  A  is  far  from  optimal, 
however,  and  it  is  hoped  that  this  might  be  improved  upon.  Interesting  open 
questions  remain  about  whether  a  finite  precision  Lanczos  computation  eventu- 
ally finds  all  eigenvalues  and  about  how  the  algorithm  behaves  if  run  for  huge 
numbers  of  steps. 

References 

[1]  T.   Chihara,    An   Introduction   to    Orthogonal  Polynomials,    Gordon   and 
Breach,  New  York,  NY,  1978. 

[2]  H.  P.  Crowder  and  P.  Wolfe,  Linear  Convergence  of  the  Conjugate  Gradient 
Method,  IBM  Jour.  Res.  and  Dev.  16:431-433,  1972. 

[3]  J.  Cullum  and  R.  Willoughby,  Lanczos  Algorithms  for  Large  Symmetric 
Eigenvalue  Computations,   Vol.  I  Theory,  Birkhauser,  Boston,  MA,  1985. 


13 


[4]  P.  Concus,  G.H.  Golub  .nd  D.P.  O'Leary,  A  Generalized  Conjugate  Gra- 
dient Method  for  the  Numerical  Solution  of  Elliptic  Partial  Differential 
Equations,  in  Sparse  Matrix  Computations,  ed.  J.R.  Bunch  and  D.J.  Rose, 
Academic  Press,  NY,  1976. 

[5]  M.  Engeli,  T.  Ginsburg,  H.  Rutishauser,  and  E.  Stiefel,  Refined  Iterative 
Methods  for  Computation  of  the  Solution  and  the  Eigenvalues  of  Self- 
Adjoint  Boundary  Value  Problems,  Birkhauser  Verlag,  Basel/Stuttgart, 
1959. 

[6]  J.  Favard,  Sur  les  Polynomes  de  Tchebicheff,  C.R.  Acad.  Sci.  Pans 
200:2052-2053,  1935. 

[7]  G.H.  Golub  and  C.F.  Van  Loan,  Matrix  Computations,  second  edition, 
Johns  Hopkins  University  Press,  Baltimore,  MD,  1989. 

[8]  J.F.  Grcar,  Analyses  of  the  Lanczos  Algorithm  and  of  the  Approximation 
Problem  in  Richardson's  Method,  Ph.D.  thesis,  University  of  Illinois,  1981. 

[9]  A.  Greenbaum,  Comparison  of  Splittings  Used  with  the  Conjugate  Gradient 
Algorithm,  Numer.  Math.  33:181-194,  1979. 

[10]  A.  Greenbaum,  Behavior  of  Slightly  Perturbed  Lanczos  and  Conjugate- 
Gradient  Recurrences,  Linear  Algebra  Appl.  113:7-63,  1989. 

[11]  M.R.  Hestenes  and  E.  Stiefel,  Methods  of  Conjugate  Gradients  for  Solving 
Linear  Systems,  J.  Res.  Nat.  Bur.  Standards  49:409-436,  1952. 

[12]  C.  Lanczos,  An  Iteration  Method  for  the  Solution  of  the  Eigenvalue  Prob- 
lem of  Linear  Differential  and  Integral  Operators,  J.  Res.  Nat.  Bur.  Stan- 
dards 45:22^-280,  1950. 

[13]  C.C.  Paige,  The  Computation  of  Eigenvalues  and  Eigenvectors  of  Very 
Large  Sparse  Matrices,  Ph.D.  thesis,  University  of  London,  1971. 

[14]  C.C.  Paige,  Error  Analysis  of  the  Lanczos  Algorithm  for  Tridiagonalizing 
asymmetric  Matrix,  J.  Inst.  Math.  Appl.  18:341-349,  1976. 

[15]  B.N.  Parlett  and  D.S.  Scott,  The  Lanczos  Algorithm  with  Selective  Or- 
thogonalization,  Math.  Comp.  33:217-238,  1979. 

[16]  J.K.  Reid,  On  the  Method  of  Conjugate  Gradients  for  the  Solution  of  Large 
Sparse  Linear  Systems,  in  Large  Sparse  Sets  of  Linear  Equations,  ed.  J.K. 
Reid,  Academic  Press,  NY,  pp.  231-254,  1971. 

[17]  H.D.  Simon,  The  Lanczos  Algorithm  with  Partial  Reorthogonalization, 
Math.  Comp.  42:115-136,  1984. 


14 


[18]  Z.  Strakos,   On  the  Real  Convergence   Rate  of  the  Conjugate  Gradient 
Method,  to  appear  in  Linear  Algebra  Appi. 

[19]  H.  Wozniakowski,  RoundofT  Error  Analysis  of  a  New  Class  of  Conjugate 
Gradient  Algorithms,  Ltn.  Alg.  Appl.  29.507-529,  1980. 


15 


IT.  -     ■     £.  :-=-  -, 

J * i * » » * «- 


:e  r"->t-.r.rr-.>  :".-«:  Differ^.:  Rho-Values 

i 1 * * i * i > * 1 i * i * 1 


1^  QBMaai  1  1  -<— < — • — • 111  I »- 


0.S  Sa'ni  I   1 — » — * * *- 


'■lII     1 *- 


s:       90      100 


103     F- 


Fig.  3.  Exact  CG  for  Different  Interval  Widths  and  F.P.  CG  for  Ax=b  (rho=.6) 


1 r 


-1 r- 


10  15 


20  25  30 

Iteration 


45 


Fig.  4.  Exact  and  Finite  Precision  CG.  n=100.  rho=.8 

10'' 

10° 

L^ 

I                            1                            I 

1                     1                    1                    1 

- 

10-3 

= 

■■.  '^^V^ 

: 

10-6 

: 

. 

- 

^\ 

If  .«.v«  ^"\^^ 

2 

10-9 

1 

= 

10-12 

- 

1 \ ; 1 

1 1 1 1 1 u 

X' 

0  10  20  30  40  50  60  70  80  90  100 

Iteration 


25 

Fig. 

5a.  Finite  Precision  Lanczos 

on  A  ( rho  = 

.6) 

X 

X 

h 

'     XX 

hi 

xxx" 

XXX 

1        ..»„        I    

xxxx          xxxx 

xxxxx 

xxxkx 

XXXXXJ 

X 

X 

XX 

XX 

XX 

XXX 

XXX           xxxx 

xxxx 

xxxx 

xxxxx 

X 

X 

X 

X 

XX 

XX 

XXX 

XXX                XXX 

xxxx 

5^ 

xxxx 

X 

X 

X 

XX 

XX 

XX 

XX                    XXX 

XXX 

xxxx 

20 

- 

31 

X 

X 

X 

X 

XXX 

XX 

XXX                XXX 

XXX 

XXX 

XXX   - 

X 

X 

X 

X 

X 

XX 

XX                  XX 

XX 

XXX 

XXX 

X 

X 

X 

X 

XX 

XX                    XX 

XX 

XX 

XXX 

X 

X 

X 

X 

X 

X 

XX                  XX 

XX 

XX 

XX 

X 

X 

X 

X 

X 

X                    XX 

XX 

XX 

XX 

15 

' 

X 

x» 

f 

X 
X 

X 
X 

X                     X 
X                     X 

XX 
X 

^ 

XX     - 
XX 

X 

I 

X 

X 

X 

X                     X 

X 

X 

]Sx 

X 

X 

X                     X 

X 

X 

X 

X 

X 

X                     X 

X 

X 

X 

10 

- 

X 

X 

X 

X 

X                   X 

X 

X 

X 

I 

X                     X 

X 

X 

X 

X 

X 

X                   X 

X 

X 

X 

X 

X                     ^ 

X 

5 

I 

X 

5| 

- 

X 

X 

X                     ^ 

x^ 

X 

X       - 

X 

n 

X 

X 

1 

X 

X 

10 


15 


20 


25 
Iteration 


30 


35 


40 


45 


50 


o« 

Fip 

.5b. 

Exact  Lanczos  on  A  (  rho  = 

.6) 

o   ' 

o 

o 

o 

<!. 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

0 

o 

o 

o 

o 

0 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

0 

o 

o 

o 

o 

o 

o 

o 

o 

15 

'■ 

o 

o 

o 
o 

o 
o 

o 
o 

- 

o 

o 

o 

o 

o 

o 

o 

10; 

r 

o 

o 
o 

o 
o 
o 
o 
o 

o 
o 
o 
o 
o 

- 

5] 

r 

o 

o 
o 

oooooo 

" 

10 


15 


20 


25 
Iteration 


30 


35 


40 


45 


50 


25 

Fig. 

5c. 

Exact  Lanczos 

on  A-hat  (  rho  =  .6  ) 

X 

X 

XX 

'     XX 

ix 

xxx' 

xxx 

XXXXX] 

X 

X 

X 

XX 

XX 

XXX 

xxx 

X 

X 

X 

XX 

XX 

XX 

xxx 

xxxx 

90 

X 

X 

X 

I" 

XX 

XX 

xxx           xxx 

xxx 

xxxx 

X 

X 

X 

XX 

XX 

XX             xxx 

xxx 

xxx 

^- 

X 

X 

X 

{x 

XX 

XX                  XX 

xxx 

xxx 

X 

X 

X 

XX 

XX                  XX 

XX 

XX 

xxx 

X 

X 

X 

X 

X 

XX                  XX 

XX 

XX 

XX 

15 

X 

X 

X 

X 

X 

X                     XX 

XX 

XX 

XX 

X 

X 

X 

X 

X 

X                     X 

XX 

XX 

XX      - 

X 

X 

X 

X                     X 

X 

XX 

XX 

1 

X 

X 

X 

X                     X 

X 

X 

Jx 

X 

X 

X                     X 

X 

X 

loi 

- 

X 

i 

X 

xx 

X 

X                     X 

X 
X 

X 

X 
X 
X 

X 
X 
X 

X 

i            I 

X 

X 

X 

X 

X 

X 

X 

5 
n 

1 

1 

X 

X 

X                     ^ 
X 

X 

X 

X 

x^    - 

X 

10 


15 


20 


25 
Iteration 


30 


35 


40 


45 


50 


25 

Fig.  6a. 

Finite  Precision  Lanczos 

on  A  ( riio  = 

•8)^ 

X   ■ 

X 

'          X 

■   X 

h 

XX                   XX 

XXX                XXX 

3bcxx 

xxx3( 

xxxx 

X 

X 

X 

XX 

XX                   XX 

XX                   XXX 

XXX 

XXX 

xxxx 

X 

X 

X 

X 

X 

XX                   XX 

XX                    XXX 

XXX 

XXX 

XXX 

X 

X 

X 

X                     XX 

XX                   XX 

XX 

XXX 

XXX 

20 

X 

X 

X 

X 

X                       XX 

XX                   XX 

XX 

^ 

XXX    - 

X 

X 

X 

h      f 

XX                 XX 

XX 

XXX 

X 

X 

X 

XX                   XX 

JSx 

XX 

XXX 

X 

X 

X 

X                     X 

X                     XX 

XX 

XX 

X 

X 

X 

X 

X                     X 

X                       XX 

XX 

XX 

XX 

15 

X 

X 

X 

X                     X 

X                     X 

XX 

XX 

XX      - 

X 

X 

X 

X                     X 

X                       X 

XX 

XX 

XX 

X 

X 

X                     X 

X                     X 

X 

XX 

XX 

X 

X 

X 

X                     X 

X                     X 

X 

]^ 

XX 

X 

X 

X                   X 

X                     X 

X 

XX 

10 

- 

X                     X 

X 

X 

X 

X 

X                     X 

X                     X 

X 

X 

XX 

X                     X 

X                     X 

X 

X 

X 

* 

X 

X 

I         5- 

X                     X 

X 

X 

X 

X 

iX                 X 
X                   X 
X                     X 

X 

X 

X 

5 

■ 

X 

»        i 

X 
X 

X 
X 

X        - 
X 

X 

X 

X 

X 

X 

X 

X                   X 

X 

X 

X 

n 

, 

-1 

' 

■ 

,         X                   X 

» 

X 

X 

10 


15 


20 


25 
Iteration 


30 


35 


40 


45 


50 


25 

Fip 

.6b. 

Exact  Lanczos  on  A  (  rho  =  .8  ) 

o    ' 

o 

o 

o 

A 

o    '                            '                            ' 

' 

0 

o 

o 

o 

0 

o 

0 

0 

o 

o 

o 

20 

■ 

o 

o 

o 
o 
o 

o 
0 
o 

o 
o 
o 

0 
0 

o 

- 

0 

o 

o 

o 

o 

o 

o 

o 

T^ 

[. 

o 

0 

o 
o 

o 
o 

o 
o 

_ 

o 

o 

0 

o 

o 

o 

o 

0 

10 

: 

o 
o 

o 

o 
o 

o 
o 
o 
o 

- 

o 

o 

o 

o 

o 

5 
n 

' 

o 

o 
o 

ooooo 

10 


15 


20 


25 
Iteration 


30 


35 


40 


45 


50 


25 

Fig. 

6c. 

Exact  Lanczos 

on  A -hat 

(rho 

==•8) 

X     ' 

X 

X 

'  t 

ix 

XX 

XXX      ' 

XXX 

■     XXX 

ijxx 

xxxSi 

xxxx 

X 

X 

XX 

XX 

XX 

XX 

XXX 

XXX 

XXX 

xxxx 

X 

X 

X 

X 

X 

XX 

XX 

XX 

XXX 

XXX 

jro? 

XXX 

20 

X 

X 

X 

X 

XX 

XX 

XX 

XXX 

XXX 

X 

X 

X 

X 

X 

XX 

XX 

XX 

XX 

XX 

XXX 

XXX    - 

X 
X 

X 

X 

X 

X 

XX 

XX 

XX 

XX 

XX 

XXX 

X 

X 

X 

X 

X 

ir 

XX 

XX 

XX 

XXX 

X 

X 

X 

X 

X 

XX 

XX 

XX 

XX 

15 

X 

X 

X 

X 

X 

X 

X 

XX 

XX 

XX 

XX 

X 

X 

X 

X 

X 

X 

X 

XX 

XX 

XX      - 

X 

X 

X 

X 

X 

X 

X 

XX 

XX 

XX 

X 

X 

X 

X 

X 

X 

X 

XX 

Sx 

X 

X 

X 

X 

X 

X 

X 

X 

f 

10 

. 

X 

X 

X 

X 

I 

X 
X 

X 
X 

XX 
X        - 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

5 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

n 

. 

, 

, 

X 

X 

} 

^  1 

X 

10 


15 


20 


25 
Iteration 


30 


35 


40 


45 


50 


m 


'><> 

Fig.  7a.  Finite  Precision  Lanczos 

on  A  f  rho  = 

1)^ 

' 

X 

X 

X 

i 

X 

XX 

XX 

XX 

ix 

XX 

XX 

X 

X 

X 

r 

XX 

XX 

XX 

XX 

X 

X 

X 

f 

XX 

XX 

XX 

X 

X 

X 

?* 

XX 

XX 

XX 

X 

X 

X 

XX 

XX 

XX      - 

X 

X 

i 

r 

X 
X 

X 
X 

XX 
X 

XX 

X 

»* 

XX 
XX 

XX 
XX 

X 

X 

X 

X 

XX 

XX 

XX 

15 

. 

X 

X 
X 

X 
X 

X 
X 

Jx 

X 
X 

X 
XX 

XX 
XX     - 

X 

X 

t 

X 

X 

X 

X 

X 

XX 

XX 

XX 

10 

X 

X 

X 

X 

X 
X 

X 
X 

X 
X 

^ 

t 

XX 

■ 

X 

X 

X 
X 

X 
X 

X 
X 

X 
X 

h 

X 
XX 

XX      - 
XX 

X 

X 

X 

X 

XX 

XX 

XX 

5 

X 

X 

X 
X 

X 
X 

X 

X 

X 
X 
X 

X 
X 
X 

XX 
X 
X 

X 

k 

XX 
XX 

X 

X 

i' 

t 

XX 

XX 

XX 

X 

X 

XX 

XX 

XX 

XX 

X 

^ 

X 

X 

^' 

XX 

XX 

XX 

XX 

n 

, 

■ 

» 

AX 

,  " 

p 

XX. 

XX 

10 


15 


20 


25 
Iteration 


30 


35 


40 


45 


50 


u 


25 

Fig.  7b. 

Exact  Lanczos  on  A  (  rho  =  1 ) 

o 

^o 

<!> 

o 

o 

o 

o 

o 

'70 

I 

o 

o 

o 

o 

0 

8 

(P 

o 

8 

0 

o 

0 

o 

15 

: 

o 

0 

0 
0 

8 

oooo 

o 

o 

o 

10 

I 

0 

o 

0 

o 

o 

o 

5 

I 

o 

o 

o 

0 

0 

0 

o 

o 

o 

o 

o 

0 

o 

o 

n 

? 

o   , 

10 


15 


20 


25 
Iteration 


30 


35 


40 


45 


50 


25 

"ig.  7c.  Exact  Lanczos 

on  A-hat  ( rho 

f  ^) 

' 

X 

X 

'     X 

i 

X 

XX       ■ 

XX 

'    XX 

ix 

XX 

XX 

X 

X 

X 

X 

XX 

XX 

XX 

XX 

XX 

XX 

X 

X 

X 

X 

XX 

XX 

XX 

XX 

XX 

20 

X 

X 

X 

X 

X 

X 

XX 

XX 

XX 

X 

X 

X 

X 

t 

i 

XX 

XX      - 

X 

X 

I 

t 

X 
X 

X 
X 

^ 

XX 
XX 

XX 
XX 

X 

X 

X 

X 

X 

X 

X 

X 

XX 

XX 

15 

. 

X 

X 

X 
X 

X 
X 

X 
X 

i^ 

X 
X 

X 

1 

XX 
XX      - 

X 

X 

i 

X 
X 

X 
X 

X 
X 

X 
X 

XX 
X 

XX 
XX 

X 

X 

X 

X 

X 

X 

X 

X 

X 

t 

f 

XX 

10 

X 

X 

X 

X 

X 

XX 

■ 

X 

X 

X 

X 
X 

X 
X 

X 
X 

h 

X 
X 

t 

XX      - 
XX 

X 

X 

X 

X 

X 

X 

XX 

XX 

X 

X 

X 

X 

X 

Jx 

XX 

XX 

XX 

5 

. 

X 

X 

X 

X 
X 

X 
X 

X 

X 

X 

X 

h 

l" 

XX 
XX 

XX 
XX      - 

X 

X 

X 

X 

X 

f 

XX 

XX 

XX 

XX 

X 

X 

h 

XX 

XX 

XX 

XX 

X 

X 

X 

X 

XX 

XX 

XX 

XX 

XX 

n 

, 

X 

} 

XX        , 

XX 

,  ^ 

p 

^. 

XX 

10 


15 


20 


25 
Iteration 


30 


35 


40 


45 


50 


Fig.  8a.  Exact  Lanczos  on  A-hat  (rho=.6,  delta=l.d-13) 
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